//      nmf.h
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#ifndef NMF_HEADER_H
#define NMF_HEADER_H

#include <boost/date_time/posix_time/posix_time.hpp> 

#define IS_MASTER_PROCESS (p_world->rank() == 0)
#define PROCESS_RANK (p_world->rank())

template<typename MATRIX>
MATRIX& nmf_bcd_W(const MATRIX& A, const MATRIX& W, const MATRIX& H, double alpha,
		  MATRIX& A_H, MATRIX& H_T_H, MATRIX& Wnew) {
  A_H = A*arma::trans(H);
  H_T_H = H*arma::trans(H);

  for (uint i = 0; i < W.n_cols; i++) {
    double h_ii = H_T_H(i,i);
    vec_t w_i = col(A_H,i)-Wnew*col(H_T_H,i);
    w_i /= (h_ii+alpha);
    w_i += col(Wnew,i)*h_ii/(h_ii+alpha);

    for (uint j = 0; j < W.n_rows; j++)
      if (w_i(j) <= 1e-16) w_i(j) = 1e-16;
    col(Wnew,i) = w_i;
  }

  return Wnew;
}

template<typename MATRIX>
MATRIX& nmf_bcd_H(const MATRIX& A, const MATRIX& W, const MATRIX& H, double alpha,
		 MATRIX& A_W, MATRIX& W_T_W, MATRIX& Hnew) {
  A_W = arma::trans(A)*W;
  W_T_W = arma::trans(W)*W;
 
  for (uint i = 0; i < H.n_rows; i++) {
    double w_ii = W_T_W(i,i);

    vec_t h_i = col(A_W,i)-arma::trans(Hnew)*col(W_T_W,i);
    h_i /= (w_ii+alpha);
    h_i += arma::trans(row(Hnew,i))*w_ii/(w_ii+alpha);

    for (uint j = 0; j < H.n_cols; j++)
      if (h_i(j) <= 1e-16) h_i(j) = 1e-16;
    row(Hnew,i) = arma::trans(h_i);
  }

  return Hnew;
}

void reset_seed();

void assign_rows(int M, std::vector<int>& begin, std::vector<int>& end_exclusive);

void read_data(matrix_t& A, int begin, int end_exclusive);

void nmf_bcd();
void master_slave_nmf();
void decentralized_nmf();

boost::posix_time::ptime now();
boost::posix_time::time_duration duration();
boost::posix_time::time_duration duration_all();
#define TOTAL_SECONDS ((double) (duration_all().total_milliseconds())/1000)
#endif
